home *** CD-ROM | disk | FTP | other *** search
/ Aminet 2 / Aminet AMIGA CDROM (1994)(Walnut Creek)[Feb 1994][W.O. 44790-1].iso / Aminet / util / misc / Fudgit233.lha / Source / src / fit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-16  |  5.3 KB  |  239 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #ifndef NOUNISTD_H
  4. #include <unistd.h>
  5. #endif
  6. #include "head.h"
  7.  
  8. #define ERRR (-1)
  9.  
  10. static double gammq(double a, double x);
  11. static double betai(double a, double b, double x);
  12. static double betacf(double a, double b, double x);
  13. static void gser(double *gamser, double a, double x, double *gln);
  14. static void gcf(double *gammcf, double a, double x, double *gln);
  15. static double sqrarg;
  16. #define SQR(a) (sqrarg=(a),sqrarg*sqrarg)
  17.  
  18. void Ft_fit(double *x, double *y, int ndata, double *sig, int mwt, double *a, double *b, double *siga, double *sigb, double *chi2, double *q)
  19. {
  20.     int i;
  21.     double wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat;
  22.  
  23.     *b=0.0;
  24.     if (mwt) {
  25.         ss=0.0;
  26.         for (i=1;i<=ndata;i++) {
  27.             wt=1.0/SQR(sig[i]);
  28.             ss += wt;
  29.             sx += x[i]*wt;
  30.             sy += y[i]*wt;
  31.         }
  32.     } else {
  33.         for (i=1;i<=ndata;i++) {
  34.             sx += x[i];
  35.             sy += y[i];
  36.         }
  37.         ss=ndata;
  38.     }
  39.     sxoss=sx/ss;
  40.     if (mwt) {
  41.         for (i=1;i<=ndata;i++) {
  42.             t=(x[i]-sxoss)/sig[i];
  43.             st2 += t*t;
  44.             *b += t*y[i]/sig[i];
  45.         }
  46.     } else {
  47.         for (i=1;i<=ndata;i++) {
  48.             t=x[i]-sxoss;
  49.             st2 += t*t;
  50.             *b += t*y[i];
  51.         }
  52.     }
  53.     *b /= st2;
  54.     *a=(sy-sx*(*b))/ss;
  55.     *siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
  56.     *sigb=sqrt(1.0/st2);
  57.     *chi2=0.0;
  58.     if (mwt == 0) {
  59.         for (i=1;i<=ndata;i++)
  60.             *chi2 += SQR(y[i]-(*a)-(*b)*x[i]);
  61.         *q=1.0;
  62.         sigdat=sqrt((*chi2)/(ndata-2));
  63.         *siga *= sigdat;
  64.         *sigb *= sigdat;
  65.     } else {
  66.         for (i=1;i<=ndata;i++)
  67.             *chi2 += SQR((y[i]-(*a)-(*b)*x[i])/sig[i]);
  68.         *q=gammq(0.5*(ndata-2),0.5*(*chi2));
  69.     }
  70. }
  71.  
  72. #undef SQR
  73.  
  74. static double gammq(double a, double x)
  75. {
  76.     double gamser,gammcf,gln;
  77.  
  78.     if (x < 0.0 || a <= 0.0) {
  79.         fputs("gammq: Invalid arguments in routine.\n", stderr);
  80.         return(ERRR);
  81.     }
  82.     if (x < (a+1.0)) {
  83.         gser(&gamser,a,x,&gln);
  84.         return 1.0-gamser;
  85.     } else {
  86.         gcf(&gammcf,a,x,&gln);
  87.         return gammcf;
  88.     }
  89. }
  90.  
  91. #define ITMAX 100
  92. #define EPS 3.0e-7
  93.  
  94. static void gser(double *gamser, double a, double x, double *gln)
  95. {
  96.     int n;
  97.     double sum,del,ap;
  98.  
  99.     *gln=lgamma(a);
  100.     if (x <= 0.0) {
  101.         if (x < 0.0) {
  102.             fputs("gser: x less < 0.\n", stderr);
  103.             return;
  104.         }
  105.         *gamser=0.0;
  106.         return;
  107.     } else {
  108.         ap=a;
  109.         del=sum=1.0/a;
  110.         for (n=1;n<=ITMAX;n++) {
  111.             ap += 1.0;
  112.             del *= x/ap;
  113.             sum += del;
  114.             if (fabs(del) < fabs(sum)*EPS) {
  115.                 *gamser=sum*exp(-x+a*log(x)-(*gln));
  116.                 return;
  117.             }
  118.         }
  119.         fputs("gser: variable a too large; ITMAX too small.\n", stderr);
  120.         return;
  121.     }
  122. }
  123.  
  124. static void gcf(double *gammcf, double a, double x, double *gln)
  125. {
  126.     int n;
  127.     double gold=0.0,g,fac=1.0,b1=1.0;
  128.     double b0=0.0,anf,ana,an,a1,a0=1.0;
  129.  
  130.     *gln=lgamma(a);
  131.     a1=x;
  132.     for (n=1;n<=ITMAX;n++) {
  133.         an=(double) n;
  134.         ana=an-a;
  135.         a0=(a1+a0*ana)*fac;
  136.         b0=(b1+b0*ana)*fac;
  137.         anf=an*fac;
  138.         a1=x*a0+anf*a1;
  139.         b1=x*b0+anf*b1;
  140.         if (a1) {
  141.             fac=1.0/a1;
  142.             g=b1*fac;
  143.             if (fabs((g-gold)/g) < EPS) {
  144.                 *gammcf=exp(-x+a*log(x)-(*gln))*g;
  145.                 return;
  146.             }
  147.             gold=g;
  148.         }
  149.     }
  150.     fputs("gcf: variable a too large; ITMAX too small.\n", stderr);
  151. }
  152.  
  153. #undef ITMAX
  154. #undef EPS
  155.  
  156. #define TINY 1.0e-20
  157.  
  158. void Ft_pearsn(double *x, double *y, int n, double *r, double *prob, double *z)
  159. {
  160.     int j;
  161.     double yt,xt,t,df;
  162.     double syy=0.0,sxy=0.0,sxx=0.0,ay=0.0,ax=0.0;
  163.  
  164.     for (j=1;j<=n;j++) {
  165.         ax += x[j];
  166.         ay += y[j];
  167.     }
  168.     ax /= n;
  169.     ay /= n;
  170.     for (j=1;j<=n;j++) {
  171.         xt=x[j]-ax;
  172.         yt=y[j]-ay;
  173.         sxx += xt*xt;
  174.         syy += yt*yt;
  175.         sxy += xt*yt;
  176.     }
  177.     *r=sxy/sqrt(sxx*syy);
  178.     *z=0.5*log((1.0+(*r)+TINY)/(1.0-(*r)+TINY));
  179.     df=n-2;
  180.     t=(*r)*sqrt(df/((1.0-(*r)+TINY)*(1.0+(*r)+TINY)));
  181.     *prob=betai(0.5*df,0.5,df/(df+t*t));
  182. }
  183.  
  184. #define ITMAX 100
  185. #define EPS 3.0e-7
  186.  
  187. static double betacf(double a, double b, double x)
  188. {
  189.     double qap,qam,qab,em,tem,d;
  190.     double bz,bm=1.0,bp,bpp;
  191.     double az=1.0,am=1.0,ap,app,aold;
  192.     int m;
  193.  
  194.     qab=a+b;
  195.     qap=a+1.0;
  196.     qam=a-1.0;
  197.     bz=1.0-qab*x/qap;
  198.     for (m=1;m<=ITMAX;m++) {
  199.         em=(double) m;
  200.         tem=em+em;
  201.         d=em*(b-em)*x/((qam+tem)*(a+tem));
  202.         ap=az+d*am;
  203.         bp=bz+d*bm;
  204.         d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
  205.         app=ap+d*az;
  206.         bpp=bp+d*bz;
  207.         aold=az;
  208.         am=ap/bpp;
  209.         bm=bp/bpp;
  210.         az=app/bpp;
  211.         bz=1.0;
  212.         if (fabs(az-aold) < (EPS*fabs(az)))
  213.             return(az);
  214.     }
  215.     fputs("BETACF: a or b too big, or ITMAX too small.\n", stderr);
  216.     Ft_catcher(ERRR);
  217.     return(ERRR);
  218. }
  219.  
  220. #undef ITMAX
  221. #undef EPS
  222.  
  223. static double betai(double a, double b, double x)
  224. {
  225.     double bt;
  226.  
  227.     if (x < 0.0 || x > 1.0) {
  228.         fputs("BETAI: Bad x.\n", stderr);
  229.         Ft_catcher(ERRR);
  230.     }
  231.     if (x == 0.0 || x == 1.0) bt=0.0;
  232.     else
  233.         bt=exp(lgamma(a+b)-lgamma(a)-lgamma(b)+a*log(x)+b*log(1.0-x));
  234.     if (x < (a+1.0)/(a+b+2.0))
  235.         return(bt*betacf(a,b,x)/a);
  236.     else
  237.         return(1.0-bt*betacf(b,a,1.0-x)/b);
  238. }
  239.